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Abstract 

The influence of periodic boundary conditions (implicit finite-size effects) 
on the anisotropy of pair correlations in computer simulations is studied for 
a dense classical fluid of pair-wise interacting krypton atoms near the triple 
point. Molecular dynamics simulation data for the pair distribution function 
gjv(r) = gN(,r,9,4>) °f -ZV-particle systems, as a function of radial distance r, 
polar angle 9, and azimuthal angle 4>, are compared directly with correspond- 
ing theoretical predictions [L. R. Pratt and S. W. Haan, J. Chem. Phys. 74, 
1864 (1981)]. For relatively small systems of N = 60, 80, and 108 atoms, 
significant angular variation is observed, which is qualitatively, and in several 
cases quantitatively, well predicted by theory. Finite-size corrections to the 
spherically-averaged radial distribution function <?jv(r), however, are found to 
be comparable to random statistical errors for runs of 10 5 time steps. 

1 Dedicated to Professor W. Gotze on the occasion of his 60 th birthday. 
PACS numbers: 61.20.Ja, 61.20.-p, 05.20.-y 
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I. INTRODUCTION 



Beyond the usual random statistical errors associated with averaging over a limited 
sample of particles in a computer simulation, systematic errors also may arise due to the finite 
size of the model system. In conventional molecular dynamics or Monte Carlo simulations of 
simple atomic systems, the size of the system usually may be chosen sufficiently large that 
finite-size effects can be safely neglected. In contrast, simulations of complex molecular (e.g., 
polymeric or amphiphilic) systems and ab initio simulations (e.g., for electronic structure) 
must often be limited to relatively small systems, for which size corrections may be a practical 
concern. 

Two general types of finite-size effect have been identified (1) explicit (or ensemble) 
size effects, caused by the suppression of density fluctuations upon fixing the number of 
particles (as in the canonical, microcanonical, and molecular dynamics ensembles); and (2) 
implicit (or anomalous) size effects, resulting from breaking of orientational symmetry due 
to imposed (usually periodic) boundary conditions. Periodic boundary conditions essentially 
eliminate surface effects in simulations of bulk phases, but entail spurious correlations, which 
can measurably affect the structure of sufficiently small systems. Both types of size effect 
are manifested in the detailed form of the pair distribution function, and thus also in any 
related thermodynamic properties. In particular, explicit size effects alter the long-range 
tail of the radial distribution function ^Ar(r) for an iV-particle system, whereas implicit size 
effects induce anisotropy in pair correlations, as reflected in angular dependence of the pair 
distribution function. 

In previous work ||, we have proposed a new method to correct for explicit finite-size 
effects in extracting the static structure factor S(Q), especially in the problematic \ow-Q 
(long-wavelength) regime, from simulation data for simple fluids. That method is based on 
the known asymptotic form of the radial distribution function gwfr) — 1 — S(0)/N, 

and on the Fourier transform relation between <7jv(r) and S(Q). Here we focus rather on 
implicit finite-size effects in simulations, and investigate their dependence on system size. 
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Some 15 years ago, a detailed formal theory (and practical approximation) for the effects of 
periodic boundary conditions on equilibrium properties of simulated fluids was put forth by 
Pratt and Haan || . At the same time, an initial, though limited, test was conducted |J using 
the available molecular dynamics (MD) data of Mandell J!]] for the anisotropic structure of a 
dense Lennard- Jones argon fluid. Otherwise, however, no further examination of the theory 
appears to have been documented. 

This motivates the present work, the main purpose of which is to present an extensive 
numerical examination of the effects of periodic boundary conditions on pair correlations in 
computer simulations. To this end, we have performed a series of standard MD simulations 
of relatively small samples of pair-wise interacting liquid krypton, a classical simple atomic 
system. Departing from standard practice, however, we have analysed the particle coor- 
dinates to extract the detailed angular dependence of the pair distribution function. For 
comparison, we have also numerically implemented the theory of Pratt and Haan ||. For 
the systems considered, significant anisotropy in pair correlations is observed, in generally 
good agreement with predictions of the theory. Although the results pertain specifically to 
atomic liquids, our conclusions are sufficiently general to be of relevance for more complex 
simulations. 

In the next section we first briefly outline the theory [||] of implicit finite-size effects, and 
then point out an implication for the scaling of such effects with system size. In Sec. Ill, 
after describing the technical details of our MD simulations and data analysis, we present 
numerical results for angular variation of the pair distribution function and compare with 
theoretical predictions. We derive as well the implied system size dependence of the more 
commonly studied radial distribution function. Finally, in Sec. IV we summarize and close 
with remarks bearing on the relative importance of implicit finite-size effects in computer 
simulations of more complex molecular systems. 
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II. THEORY OF IMPLICIT FINITE-SIZE EFFECTS 



A formally exact theory of implicit finite-size effects in computer simulations has been 
developed by Pratt and Haan (PH) |J. The theory is based on the observation that a 
periodically replicated ensemble of simulation cells is completely equivalent to an infinite 
system of oriented "supermolecules" , each consisting of a single physical (or core) atom from 
the primary simulation cell rigidly connected to all of its periodic images. Relative to the 
core atom, the image atoms sit on the sites of a lattice whose unit cell is the simulation cell. 
The key point is that equilibrium statistical mechanical properties of the finite simulated 
system may be described by precisely the same theoretical methods commonly applied to 
bulk molecular systems. In particular, the pair distribution function #Ar(r) for the iV-particle 
system can be related to its equivalent for the infinite supermolecular system. For an atomic 
system with uniform one-particle density p(r) = p, the pair distribution function is defined 
such that {pg^{r)dr) is equal to the ensemble average of the number of atoms in a volume 
element df at position r, in the presence of a particle fixed at the origin. Thus, gjv(r) is 
proportional to the conditional probability of finding a particle at r, given a particle at the 
origin. Using standard cluster expansion techniques |8]||, In ^(r) can be formally expressed 
as an exact expansion in terms of two-supermolecule Mayer cluster functions ||. 

If interactions are of sufficiently short range that a given particle does not interact with 
any of its periodic images, as is the case when the pair potential is truncated within the 
simulation cell, then an important class of "spring bond" graphs can be summed to all 
orders f|. Working in the grand canonical ensemble to avoid explicit size effects, and 
summing exactly a certain class of "unbridged" graphs, which depend on only the infinite- 
system function g(r), but neglecting another class of "bridged" graphs, PH have derived 
a practical approximate relation between gAr(f) and g(r), which may be expressed in the 
superposition-like form 

9 N(r 12) - g(r 12 )l[g{\fi -r 2i |), (1) 
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where r\ 2 = r i — f 2 and the product index runs over all periodic images of particle 2. 
Physically, Eq. (|1|) fulfills the intuitive expectation that the statistical correlation of a given 
particle 1 with a second particle 2 depends not only on the position of particle 2, but also on 
the positions of all periodic images of particle 2. Evidently, gfjv(f) is anisotropic to the extent 
that distances between image particles fall appreciably within the range of pair correlations 
in the corresponding infinite system. Expansion ([!]) is exact to first order in the bulk 
density and, since the neglected graphs are relatively highly connected, the approximation 
is expected to be accurate for sufficiently large system sizes and short-ranged interactions. 
Numerical computation of gN{r) from Eq. (|l|) clearly requires prior knowledge of the bulk 
function g(r). Practical implementation of the theory is discussed in Sec. Ill B. 

The characteristic dependence of gN(r) on system size N implied by the PH theory can 
be inferred by the following argument. The size dependence clearly arises because the image 
distances \f\ — r 2i \ depend implicitly on the cell length L. Suppose that L is large enough 
that the shortest image distance already lies in the long-range asymptotic tail of g(r). Then, 
in terms of the pair correlation function h(r) = g(r) — 1, we can write [from Eq. ([]])], 



g N (r 12) /g(r 12 ) = n 



1 + h(\n - f 2i \) 



~l + J2h(\n-r 2i \), (2) 

i 

to leading order in h{r). Therefore, the system size dependence of gjv(r 12 ) is dictated by the 
asymptotic behaviour of h(r). Suppose, for example, that h(r) has power-law asymptotic 



behaviour pi of the form h(r) ~ l/r n (r — > 00). Since each term in the sum then scales 



with system size as 1/L n (to leading order), and since L ~ iV 1//3 , we conclude that 

gN(r 12 )/g(r l2 ) ~ l + 0(7-M n/3 V (3) 



Thus, assuming power-law asymptotic pair correlations, the PH theory predicts that for 
n > 3 implicit finite-size corrections to g7v(r) decrease faster with system size than do 
explicit 0(1 /N) corrections. Of course, for more rapidly (e.g., exponentially) decaying pair 
correlations the size corrections would decrease correspondingly faster. 



III. MOLECULAR DYNAMICS SIMULATIONS 



A. SIMULATIONS AND ANALYSIS 

In order to test the theory described in Sec. II, we have performed isothermal-isochoric 
MD simulations for systems of N = 60, 80, and 108 atoms of krypton, a well-studied 



experimental system ||TT|. Periodic boundary conditions were applied to a cubic simulation 
cell of length L, such that a particle leaving the cell through one face was simultaneously 
replaced by an image particle entering through the opposite face. The atoms were assumed 



to interact via the Aziz pair potential |T2|, which is known to give an accurate description 
of the structure and thermodynamics of rare-gas systems. For krypton, the pair potential 
is parametrized by the distance at which it crosses zero, o = 3.579A, the distance at 
which it attains a minimum, a m = 4.012A, and the well depth e/ks = 200K. To avoid 
direct interaction between a given particle and its periodic images, the pair potential was 
truncated and shifted to zero at a cut-off distance r c = L/2. Particle trajectories were 
computed by solving the classical equations of motion with a fifth-order Gear's predictor- 
corrector algorithm. The thermodynamic state for all of our simulations is defined by the 
reduced density p* = pal = 0.8 and the reduced average temperature T* = ksT/e = 0.7, 
placing the system near its triple point. Aside from an initial equilibration phase, during 
which velocities were rescaled to establish the desired average temperature, the total number 
of particles, volume, energy, and momentum were kept fixed. Runs of 10 5 time steps, with a 
time step of 0.005 ps, generated about 1000 independent configurations, sufficient to ensure 
that random statistical errors are small in comparison with the finite-size effects of interest. 

As a preliminary test of the simulation and analysis methods, we first compute the 
radial distribution function (?7v(r), defined such that the quantity pg7v(r)(47rr 2 Ar) is equal 
to the ensemble average of the number of atoms in a spherical shell of thickness Ar at radial 
distance r away from a particle centred at the origin (r = 0). Note that the isotropic function 
^Ar(r) is the spherical average of the anisotropic function gN{r) defined in Sec. II. Following 
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standard practice, QNij) is computed by binning particles in radial shells centred on a given 
central particle and then averaging over all central particles and over all configurations. 

Figure 1 shows the resulting ^jv(^) for the system size N = 108 out to a distance of r = 
L/2. Choosing the shell thickness Ar = 0.028cx o , the random statistical errors are at most 
±0.02. Corresponding results for N = 60 and N = 80 are found to be practically identical 
to within statistical error, suggesting that the function g^fr) is a reasonable measure of the 
bulk g(r), and that at this level finite-size effects are relatively minor. We return to this 
point below. Also shown in Fig. 1, for comparison, is the prediction of the perturbation 
theory of Weeks, Chandler, and Andersen (WCA) Jl3| , according to which 



g(r) ~ exp[-u (r)/k B T]y PY (r/d;pd 3 ). (4) 

Here u a (r) is the repulsive part of the pair potential, ypy is the Percus-Yevick (PY) indi- 
rect correlation (or cavity) function of the hard-sphere fluid, and d is the WCA prescrip- 
tion for the effective hard-sphere diameter. At the thermodynamic state investigated, 
d = 1.027cx o . For ypy, we use the analytic solution of the PY equation for hard spheres P,|H| . 
The WCA prediction is seen to be in reasonable agreement with the simulation data, despite 
minor discrepancies around the first and second peaks. 

As discussed in Sees. I and II, pair correlations in a finite system are not necessarily 
isotropic [i.e., g^if) ^ gN{ r )] because of broken orientational symmetry. As a measure of 
anisotropy, we have examined the angle- dependent pair distribution functions 

1 r 1 

g N (r, 4>) = - J d(cos9)g N (r) (5) 

and 

1 r 27T 

g N (r,9) = — d<j)g N (r), (6) 

271 JO 

where pjv(r) = <7jv(V, 9, (ft), with r, 9, and (ft being the radial distance, polar angle, and 
azimuthal angle. Note that the radial distribution function is obtained by integrating either 
Eq. (H) with respect to or Eq. @ with respect to cos 9. By convention, the Cartesian 



axes specified by {9 = tt/2, = 0), (9 = n/2, <fi = vr/2), and (9 = 0,(f) arbitrary) are directed 
perpendicular to the sides of the simulation cell. Symmetry dictates that <?jv(r, 0) is periodic 
in with period tt/2 and is mirror symmetric about = tt/A. Similarly, g7v(r, 9) is periodic 
in 9 with period 7r and is mirror symmetric about 9 = tt/2. All angles therefore can be 
represented in the reduced ranges < < 7r/4 and < 9 < tt/2. In practice, gjy(r, 0) is 
computed by counting the number of particles within a particle- centered longitudinal wedge 
of longitude 0, radius r, angular width A0, and radial thickness Ar, and then averaging 
over all central particles and all configurations. Similarly, gjy(r,9) is computed by counting 
the number of particles within a latitudinal slice of latitude 9, radius r sin 9, angular width 
A9, and radial thickness Ar, averaged over all central particles and configurations. 



B. RESULTS AND COMPARISON WITH THEORY 

Figures 2-4 display our MD results for the angle-dependent distribution function gjv(r, 0), 
as a function of azimuthal angle [Eq. (|^)], at fixed radial distances rjo Q = 1.1, 1.6, and 2.1 - 
near the first maximum, the first minimum, and the second maximum of <?Ar(r), respectively 
(see Fig. 1). For the smallest system (Fig. 4), only the data for the first two distances are 
shown, since the cut-off at r = L/2 occurs before the second maximum of g7v(r). Figure 5 
presents similar results (largest system only) for g^(r,9) as a function of polar angle 9 
[Eq. (§)]. Note that as 9 — > the statistical error increases because of the geometrical factor 
sin9 in the volume of the latitudinal slice. In order to accumulate reasonable statistics, 
we have used angular widths A0 = A9 = 7r/36 and a radial thickness Ar = 0.112cx o . As 
the latter is four times larger than the shell thickness used in the calculation of ^^(r), the 
average values of (^(r, 0) and <?jv(r, 9) differ slightly from the corresponding values of <?jv(r). 
The angular variations illustrated in Figs. 2-5 are clearly significant - up to five times larger 
than the statistical error in the case of 0, and up to ten times larger in the case of 9. We 
mention in passing that we have also analysed the one-particle density p(r) and confirmed 
it to be translationally uniform, i.e., p(r) = p. 



S 



Also shown in Figs. 2-5, for comparison, are corresponding predictions of the PH theory. 
As noted in Sec. II, numerical implementation of the theory requires knowledge of the bulk 
system radial distribution function g(r). To allow for a consistent comparison between 
theory and simulation, for distances r < L/2 ~ 2.56a Q we approximate the required infinite- 
system g(r) by its finite-system counterpart gN{ r ) computed from our N = 108 simulation. 
This represents a good approximation to the true bulk function, as shown in Sec. II and 
discussed further below. For larger distances r > L/2, beyond the range of the simulation 
data, we use the g{r) predicted by the WCA perturbation theory [|T3|1 . As demonstrated in 



Sec. II, this gives a reasonable fit for the relatively dense system considered here (see Fig. 1). 
(We note, however, that at lower densities, where perturbation theories tend not to perform 
as well, more accurate integral-equation theories may be required.) The function <?jv(r) is 
computed directly from Eq. ([!]), taking into account all image particles in a 5x5x5 ensemble 
of cells centred on the primary cell (124 images). Numerical integration with respect to 
9 or cf) then yields the distribution functions defined by Eqs. (0) and (|6]). Evidently, the 
PH theory tracks the angular variations remarkably well, at least at the thermodynamic 
state that we have simulated. Where quantitative discrepancies occur, the theory tends to 
underestimate the magnitude of the variations. As a check on the calculations, note that 
(?jv(r, = 0) = <7iv(r, 6 = vr/2), as should be so by symmetry. 

Finally, we examine the effect of periodic boundary conditions on the radial distribution 
function by computing the PH prediction for the spherically-averaged function 

1 r 2n 



1 f f 

8gN(r) = — / d(f) d(cos6)g N (r) 
Atc Jo J-i 



9(r). (7) 



This gives a measure of the deviation of the iV-particle function ^Ar(r) from its bulk coun- 
terpart g(r). For simplicity, we use in this case the WCA g(r) for all distances. As Fig. 6a 
illustrates, the predicted 5giy(r) oscillates as a function of r about zero with an amplitude 
of less than 0.02 for iV = 108, increasing as the system size decreases. The oscillations arise 
from the fact that as r increases the image distances |r*i — r 2 j| in Eq. ([I]) alternate between 
regions where g(r) is greater or less than unity. The amplitude is in fact comparable to the 



random statistical errors in our simulations (see Fig. 1), although of course for longer runs it 
would begin to exceed statistical error. Interestingly, at distances near the first maximum, 
first minimum, and second maximum [rja = 1.1, 1.6, and 2.1) the deviation from bulk 
behaviour essentially vanishes. This further justifies our use above of the MD ^Ar(r) as an 
approximation to the bulk g(r) in numerically implementing the PH theory to study angular 
variations. It should be emphasized that the practical utility of the PH theory for predicting 
finite-size effects on thermodynamic properties related to volume integrals of g(r), such as 
the equation of state, is limited by an inherent thermodynamic inconsistency in the approx- 
imation of Eq. ([!]) • In this connection, however, exact implicit finite-size corrections to 



virial coefficients have been derived fl6|| . 

The predicted scaling of (?at(V) with system size is shown in Fig. 6b at two fixed distances. 
The oscillations arise from the alternation of image distances, with increasing L, between 
regions where g(r) > 1 and regions where g(r) < 1, and have a wavelength increasing as 
iV 1 / 3 . Clearly the iV-particle function decays rapidly to its bulk limit with increasing N. 
The functional form is revealed by the log-log plot (inset) of the maxima and minima of 
the bounding envelope, indicating a near-linear relation. For the short-ranged pair potential 
used here, this implies a simple power law of the form gN{r)/g{r) ~ 1 + 0{1/N V ) (with 
v ~ 2.2), consistent with the scaling argument presented in Sec. II. 



IV. SUMMARY AND CONCLUSIONS 

In summary, we have studied the magnitude of finite-size effects induced by periodic 
boundary conditions through a series of molecular dynamics simulations of a dense krypton 
fluid of between N = 60 and N = 108 atoms near its triple point. As a measure of 
anisotropy, we have analysed the variation with azimuthal and polar angles of the pair 
distribution function at several radial distances and system sizes. Significant anisotropy 
was observed, the magnitude of angular variations being several times larger than random 
statistical errors for runs of 10 5 time steps. 
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For comparison, the theory of Pratt and Haan || was implemented and was found to sat- 
isfactorily model the anisotropies observed in the simulations. The same theory predicts that 
spherical averaging of the pair distribution function considerably diminishes the magnitude 
of implicit size corrections, reducing them to the level of our statistical errors. The theory 
further implies that for asymptotic pair correlations of the power-law form (~ l/r n ) the size 
correction scales with system size as 0(1/ N n ^ 3 ). This is confirmed by our comparison of 
simulation data for the radial distribution function at different system sizes, which indicates 
that implicit corrections decrease rapidly with system size, significantly more rapidly than 
explicit corrections. 

We conclude by pointing out that since implicit finite-size effects can obviously cause 
significant anisotropy in pair correlations even in the case of a simple atomic liquid with 
spherically symmetric pair potential, they are certainly also a potential source of error in 
extracting orientational correlation parameters from simulations of molecular systems, es- 
pecially for anisometric molecules JIT). In MD simulations of 256 CS2 molecules ||18|| , for 



instance, significant anisotropy in the Kirkwood #2 parameter was observed [|17| . Further- 
more, the distortion of three-particle and higher-order correlations by boundary conditions 
is expected to be even more pronounced || . Finally, we emphasize that here we have consid- 
ered only the case of short-range interactions. The case of long-range interactions is relevant 



to simulations of molten salts, charge-stabilized colloids, and other charged systems ||19|| . Fu- 
ture study of these issues, by simulation and by practical extension of the theory of finite-size 
effects would be worthwhile. 
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FIGURES 

FIG. 1. MD simulation data (dots) and WCA prediction |l3| (curve) for the radial distribution 
function of 108 krypton atoms at reduced density p* = 0.8 and reduced temperature T* = 0.7. 
Error bars in the inset represent random statistical errors of plus or minus one standard deviation. 

FIG. 2. MD data (symbols) compared with theoretical predictions || (curve) for the 
angle-dependent pair distribution function g]\f(r,(p) vs. azimuthal angle <j) (in radians, with 
< (j) < vr/4) for N = 108 and radial distances (a) r/a Q = 1.1, (b) 1.6, and (c) 2.1. Error 
bars represent random statistical errors of plus or minus one standard deviation. 

FIG. 3. Same as Fig. 2, but for N = 80. 

FIG. 4. Same as Fig. 2, but for N = 60. 

FIG. 5. MD data (symbols) compared with theoretical predictions [|| (curve) for the an- 
gle-dependent pair distribution function g;\i(r,9) vs. polar angle 9 (in radians, with < 9 < tt/2) 
for N = 108 and radial distances (a) r/a = 1.1, (b) 1.6, and (c) 2.1. 

FIG. 6. Theoretical predictions H for the iV-particle radial distribution function gpf(r): (a) 
deviation from the bulk g{r) vs. radial distance r < L/2 for N = 108 (solid curve), N = 80 
(long-dashed), and N = 60 (short-dashed); (b) scaling with system size N for r = 1.6cr (solid) 
and r = 2.1a (dashed); inset shows scaling of bounding envelope on a log-log scale for r = 1.6ct 
(circles) and r = 2.1cr (squares). 
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Fig. 1 Inset (Denton and Egelstaff) 




Fig. 1 (Denton and Egelstaff) 




Fig. 2b (Denton and Egelstaff) 
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Fig. 3a (Denton and Egelstaff) 




Fig. 3b (Denton and Egelstaff) 




Fig. 4a (Denton and Egelstaff) 




Fig. 4b (Denton and Egelstaff) 




Fig. 5a (Denton and Egelstaff) 




Fig. 5b (Denton and Egelstaff) 




Fig. 6b Inset (Denton and Egelstaff) 
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